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Abstract 

The flutter characteristics of the first AGARD standard 
aeroelastic configuration for dynamic response. Wing 445.6, 
are studied using an unsteady Navier-Stokes algorithm in 
order to investigate a previously noted discrepancy between 
Euler flutter characteristics and the experimental data. The 
algorithm, which is a three-dimensional, implicit, upwind 
Euler/Navier-Stokes code (CFL3D Version 2.1), was previ- 
ously modified for the time-marching, aeroelastic analysis 
of wings using the unsteady Euler equations. These modifi- 
cations include the incorporation of a deforming mesh algo- 
rithm and the addition of the structural equations of motion 
for their simultaneous time integration with the governing 
flow equations. In this paper, the aeroelastic method is ex- 
tended and evaluated for applications that use the Navier- 
Stokes aerodynamics. The paper presents a brief descrip- 
tion of the aeroelastic method and presents unsteady cal- 
culations which verify this method for Navier-Stokes cal- 
culations. A linear stability analysis and a time-marching 
aeroelastic analysis are used to determine the flutter char- 
acteristics of the isolated 45° swept-back wing. Effects of 
fluid viscosity, structural damping, and number of modes in 
the structural model are investigated. For the linear stabil- 
ity analysis, the unsteady generalized aerodynamic forces of 
the wing are computed for a range of reduced frequencies 
using the pulse transfer-function approach. The flutter char- 
acteristics of the wing are determined using these unsteady 
generalized aerodynamic forces in a traditional V-g analysis. 
This stability analysis is used to determine the flutter char- 
acteristics of the wing at free-stream Mach numbers of 0.96 
and 1.141 using the generalized aerodynamic forces gener- 
ated by solving the Euler equations and the Navier-Stokes 
equations. Time-marching aeroelastic calculations are per- 
formed at a free-stream Mach number of 1.141 using the 
Euler and Navier-Stokes equations to compare with the lin- 
ear V-g flutter analysis method. The V-g analysis, which is 
used in conjunction with the time-marching analysis, indi- 
cates that the fluid viscosity has a significant effect on the 
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supersonic flutter boundary for this wing while the structural 
damping and number of modes in the structural model have 
a lesser effect. 

Nomenclature 


Aij generalized aerodynamic force resulting from 

pressure induced by mode j acting through 
mode i 

b root semichord 

c root chordlength 

C P pressure coefficient 

(ji structural damping for mode i 

k reduced frequency, -rfp— 

M . free-stream Mach number 

Q free-stream dynamic pressure 

Re c free-stream Reynolds number based on root 

chordlength 

T dimensional time 

I ’/ ■ flutter speed 

Leo streamwise free-stream speed 

a angle of attack 

i j nondimensional semispan location 

p mass ratio 

uj angular frequency 

uj a uncoupled natural frequency of the wing first 

torsion mode 

Introduction 

From the first calculations of Ballhaus and Goorjian 1 
to more recent applications of the Euler and Navier-Stokes 
equations for the aeroelastic analysis of three-dimensional 
wings, the field of computational aeroelasticity has pro- 
gressed rapidly. 2 With the advent of more powerful com- 
puters and more efficient algorithms, researchers have in 
recent years computed more detailed and sophisticated sim- 
ulations of aeroelastic phenomena. Because of their rel- 
ative computational efficiency, the transonic small distur- 
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bance (TSD) equation and the full potential (FP) equation 
have been applied to a wider variety of three-dimensional 
configurations than the Euler and Naiver-Stokes equations. 2 
Similarly, the TSD and FP equations also have been utilized 
for more detailed analyses of the aeroelastic characteristics 
of these configurations. In contrast, the higher-order meth- 
ods, such as the Euler and Navier-Stokes equations, typi- 
cally have been applied to the analyses of flexible wings 
over a limited range of flow conditions. 3-11 Only recently, 
the Euler equations have been used to compute the flutter 
boundary for the AGARD standard aeroelastic configura- 
tion for dynamic response. Wing 445.6. 12- 13 In Ref. 12, an 
unstructured-grid Euler code is used to compute the flutter 
characteristics of Wing 445.6 for free-stream Mach num- 
bers ranging from 0.499 to 0.96, and the computed flutter 
speeds and frequencies are compared with the experimental 
values measured for this wing. In Ref. 13, Lee-Rausch and 
Batina compute the complete flutter boundary for the same 
wing using the Euler equations on a structured grid. Calcu- 
lated flutter results are compared with experimental data for 
seven free-stream Mach numbers, which define the flutter 
boundary over a range of Mach number from 0.499 to 1.14 
(See Fig. 1.) These comparisons show good agreement in 
flutter characteristics for free-stream Mach numbers below 
one. However, for free-stream Mach numbers above one, the 
computed aeroelastic results predict a premature rise in the 
flutter boundary as compared with the experimental bound- 
ary. The purpose of this paper is to extend the capability 
presented in Ref. 13 and to compute the flutter character- 
istics of the AGARD Wing 445.6 using the Navier-Stokes 
equations so that the source of the discrepancy in the flutter 
characteristics at the supersonic free-stream Mach numbers 
can be investigated. 

Many of the aeroelastic analyses performed with the 
nonlinear flow equations have been obtained by calculating 
the transient response of the coupled aerodynamic/structural 
system. Other analyses have been performed using harmonic 
loads. 14-17 In the harmonic analysis, the unsteady aerody- 
namic forces are assumed to be locally linear and a tradi- 
tional stability analysis is used to determine flutter charac- 
teristics based on the calculated harmonic loads. One advan- 
tage of the linear stability analysis is that it offers a signifi- 
cant computational savings over the time-marching analysis 
if the number of modes necessary to model the structure is 
small and if there are no static aeroelastic deformations. An- 
other advantage of the linear stability analysis is that it can 
provide more information regarding the effects of individ- 
ual modes and the effects of structural damping from fewer 
computations than the time-marching approach. It is impor- 
tant to note, however, that the assumption of superposition 
of airloads may not be accurate in some cases such as those 
involving large amplitudes of motion and separated flow. 
However, in previous studies, 14 - 15 - 17 for cases with low an- 
gles of attack and small amplitudes of motion, the flutter 
characteristics obtained from time-marching analyses com- 
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Figure 1 Comparison of computed Euler flutter results 
from Ref. 13 with experimental data for Wing 445.6. 


pare favorably with those obtained from a harmonic analy- 
sis. For these reasons, the linear stability analysis is used in 
conjunction with time-marching analyses to study the aeroe- 
lastic characteristics of the AGARD Wing 445.6 using the 
Navier-Stokes equations. 

In Ref. 13, modifications are made to an existing three- 
dimensional, unsteady Euler/Navier-Stokes code (CFL3D 
Version 2.1) for the aeroelastic analysis of wings. These 
modifications include the incorporation of a deforming mesh 
algorithm and the addition of the structural equations of mo- 
tion for their simultaneous time integration with the govern- 
ing flow equations. This paper gives a brief description 
of these modifications and presents unsteady calculations 
which verify the modifications to the code for the solution of 
the Navier-Stokes equations. Results from calculations per- 
formed about a rigid wing undergoing forced plunging and 
pitching motions are presented to verify the performance 
of the deforming mesh algorithm. A linear stability anal- 
ysis and a time-marching aeroelastic analysis are used to 
determine the flutter characteristics of the wing. Effects of 
fluid viscosity, structural damping, and number of modes 
in the structural model are investigated. For the linear sta- 
bility analysis, the unsteady generalized aerodynamic forces 
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(GAF’s) of the wing are computed for a range of reduced 
frequency using a pulse transfer-function analysis. The flut- 
ter characteristics of the wing are determined using these 
unsteady GAF’s in a traditional V-g analysis. This stability 
analysis is used to determine the flutter characteristics of the 
wing at free-stream Mach numbers of 0.96 and 1.141 using 
the GAF’s generated with the Euler equations and with the 
Navier-Stokes equations. In addition, time-marching aeroe- 
lastic calculations are performed at a free-stream Mach num- 
ber of 1.141 using the Euler and Navier-Stokes equations 
for comparison with the linear V-g flutter analysis results. 
The time-marching Euler flutter characteristics are recom- 
puted on a mesh that is similar to the Navier-Stokes mesh 
in order to more effectively isolate the effects of viscosity. 
Computed flutter results are compared with the experimental 
data and with other computational results from Ref. 13 to 
demonstrate the aeroelastic capability for the Navier-Stokes 
equations and to evaluate the effects of viscosity on the flut- 
ter characteristics at these flow conditions. 

Time-Marching Aeroelastic Analysis 

The time-marching aeroelastic procedure used in this 
study is typical of those currently in use. 11- 12 - 18 In general, 
the aeroelastic equations of motion are formulated in terms 
of a finite modal series of free-vibration modes. These equa- 
tions then are written in terms of a linear state-space equation 
such that a modified state-transition-matrix integrator can be 
used to march the coupled fluid-structural system forward in 
time. The fluid forces are coupled with the structural equa- 
tions of motion through the generalized aerodynamic forces. 
To determine the flutter conditions at a given free-stream 
Mach number, aeroelastic transients are computed at several 
values of dynamic pressure which bracket the flutter point. 
The frequency and damping characteristics of the transient 
responses at each dynamic pressure are determined from a 
least squares curve fit, 19 and the flutter dynamic pressure 
and frequency associated with this Mach number can be es- 
timated by interpolation. 

The time-marching aeroelastic procedure for a problem 
utilizing the Navier-Stokes equations is almost identical to 
the procedure described in more detail in Ref. 13 for the 
Euler equations. However, when using the Navier-Stokes 
equations, a free-stream Reynolds number must be speci- 
fied for the static and dynamic computations at each dy- 
namic pressure of interest. This requirement raises the issue 
of how to vary the dynamic pressure in a computational 
fluid dynamics, time-marching aeroelastic analysis to obtain 
the flutter conditions for a specific free-stream Mach num- 
ber. When comparing with experimental data, either the 
free-stream density or the free-stream velocity can be devi- 
ated from the experimental value to obtain variation in the 
dynamic pressure. Changing the density or velocity will, 
however, change the Reynolds number associated with that 
flow condition. An analysis of each of these options for 
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(a) Moo = 0.96. 



0 1 2 
Q/Qexp 

(b) Moo = 1.141. 

Figure 2 Variation in free-stream Reynolds 
number with variation in dynamic pressure 
for time-marching aeroelastic calculations. 

the experimental flutter conditions of the Wing 445.6 at 
Moo = 0.96 and 1.141 is shown in Figs. 2(a) and (b), re- 
spectively. Figures 2(a) and (b) show the variation in the 
free-stream Reynolds number (Re c ) with the variation in the 
free-stream dynamic pressure nondimensionalized by the ex- 
perimental flutter dynamic pressure OV/ 1 / • )• Figures 2(a) 
and (b) indicate that for a constant free-stream density, vary- 
ing the free-stream velocity results in a smaller variation in 
Reynolds number than maintaining a constant velocity while 
varying the density (especially for conditions above the ex- 
perimental dynamic pressures). Varying the free-stream ve- 
locity for a given Mach number does, however, require a 
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change in the free-stream temperature and, subsequently, the 
adiabatic wall temperature. The effect of Reynolds number 
variation on the time-marching cases considered in this study 
is discussed later in the paper. 

Linear Stability Analysis 


A conventional V-g method is used in this study for 
computing flutter speeds and frequencies. This method in- 
terpolates the given generalized aerodynamic forces (GAF’s) 
to compute the eigenvalues of the flutter determinant. The 
procedure used to obtain the flutter eigenvalues is contained 
in the parametric flutter analysis program FLUTDET which 
is part of the FAST flutter analysis package. 20 The capabili- 
ties and techniques used in FLUTDET are described in Ref. 
21. Although FLUTDET originally was designed to inter- 
face with GENFLU, which calculates unsteady aerodynamic 
forces from subsonic kernel matrices, compatible input from 
other aerodynamic methods can be used. In addition to the 
GAF’s, the generalized mass and modal stiffnesses are re- 
quired input to FLUTDET. The program generates flutter 
characteristics and V-g plots for each specified free-stream 
density. 

Traditional flutter analysis methods, such as the V-g 
method, assume that the unsteady aerodynamic forces are 
locally linear and utilize modal superposition of harmonic 
loads. Linear aerodynamic methods typically compute the 
GAF’s as functions of reduced frequency. For nonlinear 
aerodynamic methods, these GAF’s are obtained from time- 
marching calculations by computing several cycles of a 
forced harmonic oscillation and using the last cycle of os- 
cillation to determine the first harmonic component of the 
GAF’s. Multiple time-marching computations are required 
at various reduced frequencies to generate the GAF’s for 
each structural mode used in the flutter analysis. 

An alternate method for determining the GAF’s is the 
pulse transfer-function analysis. In contrast to the forced 
harmonic method, the pulse analysis can determine the 
GAF’s over a range of reduced frequency in a single time- 
marching calculation for each mode. In the pulse analy- 
sis, the unsteady forces are determined from the response 
due to motion represented by a smoothly varying, exponen- 
tially shaped pulse. A fast Fourier transform of the unsteady 
force is divided by the fast Fourier transform of the dis- 
placement to obtain the GAF. Transonic small disturbance 
results computed using the pulse analysis for a pitching flat 
plate are found to be in good agreement with linear theory 
calculations. 22 Also in Ref. 22, the GAF’s of airfoils at tran- 
sonic speeds computed from a pulse analysis are shown to be 
in good agreement with the GAF’s computed using the har- 
monic method which tends to verify that the analysis is valid 
for predicting the linear small perturbation response about 
a nonlinear flowfield. Similarly in Refs. 11 and 23, Euler 
results for a pitching and plunging airfoil show good agree- 
ment between the pulse and harmonic methods. Therefore, 


because of its computational efficiency, the pulse transfer- 
function method is used in this study to compute the GAF’s 
for input to the V-g analysis. 

For linear methods, the GAF’s for a vibration mode 
are functions of free-stream Mach number, planform (ge- 
ometry), and reduced frequency. For nonlinear methods 
such as the TSD and Euler equations, if the assumptions 
of local linearity are maintained, then this functionality also 
will hold true. However, for the Navier-Stokes equations 
an additional parameter must be considered: the free-stream 
Reynolds number. The V-g stability analysis does not ensure 
that the Reynolds number associated with the computation of 
the GAF’s will match the Reynolds number associated with 
the computed flutter condition. The only way to ensure a 
matched Reynolds number solution is to iterate between the 
calculations of the GAF’s and the V-g analysis. This option 
is not feasible for a Navier-Stokes analysis due to constraints 
on computational resources. However, for the cases ana- 
lyzed in this study, it was not considered to be a necessary 
requirement for the following reason. In the V-g analysis, 
the flutter speed is computed for a given free-stream den- 
sity. Typically when comparing with experimental data, the 
experimental free-stream density is specified. As discussed 
in the previous section, variations in free-stream Reynolds 
number due to changes in free-stream velocity are small. 
If the GAF’s are computed at the experimental Reynolds 
number, the difference between the experimental Reynolds 
number and the Reynolds number based on the computed 
flutter characteristics are small for the cases considered in 
this study. These differences are discussed later in the paper. 

Upwind Euler/Navier-Stokes Algorithm 

The CFL3D 24 20 code uses a three-factor, implicit, 
finite-volume algorithm based on upwind-biased spatial dif- 
ferencing to solve the time-dependent Euler equations and 
thin-layer Navier-Stokes equations. The algorithm, which is 
a cell-centered scheme, uses upwind differencing based on 
either flux-vector splitting or flux-difference splitting. Both 
types of upwind differencing account for the local wave- 
propagation characteristics of the flow and sharply capture 
shock waves. Also, because these schemes are naturally dis- 
sipative, additional artificial dissipation terms are not nec- 
essary. Several types of flux limiting are available within 
the code to prevent oscillations in the solution near shock 
waves which are typically found in higher-order schemes. 
For applications utilizing the thin-layer Navier-Stokes equa- 
tions, two turbulence models are available: the equilibrium, 
algebraic, eddy-viscosity model of Balwin-Lomax, 27 and the 
nonequilibrium half-equation model of Johnson and King. 28 
For unsteady cases, the original algorithm contains the nec- 
essary metric terms for a rigidly translating and rotating 
mesh which moves without deforming. For cases involving 
a deforming mesh, however, an additional term accounting 
for the change in cell volume must be included in the time- 
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discretization of the governing equations. This modification 
is implemented as described in Ref. 11. 

Deforming Mesh Algorithm 


In the time-marching aeroelastic calculations and in the 
modal pulse transfer-function calculations, the mesh must 
be updated at every time level so that it conforms to the 
aeroelastically deformed shape of the wing. Because the 
aeroelastic motions of the wing are arbitrary in nature, a 
general mesh updating procedure is necessary. One such 
method, the deforming mesh algorithm, models the mesh 
as a network of springs and solves the static equilibrium 
equations for this network to determine the new locations 
of the mesh grid points. This algorithm was originally 
developed by Batina 29 for tetrahedral cells and extended 
by Robinson et al. 11 for hexahedral cells. The method 
described in Ref. 11 was used by the current authors 
for time-marching aeroelastic calculations on Euler grids. 13 
This algorithm is extended in the current study for use on 
Navier-Stokes meshes. The basic principles of the method 
do not change although some modifications to the deforming 
mesh boundary conditions are required as described below. 

For the deforming mesh algorithm, the edge of each 
hexahedral cell is modeled as a spring with a stiffness that 
is inversely proportional to a power of the length of the edge. 
Diagonal springs are added along the faces of each cell in or- 
der to prevent cell shearing. Similarly, the stiffness of these 
springs is inversely proportional to a power of the length 
of the diagonal. As suggested in Ref. 11, a power of three 
was used in the present calculations. At each time level, the 
grid points on the outer boundary are held fixed, and the 
displacement of the wing surface is specified. For aeroe- 
lastic calculations, the displacement is determined from the 
integration of the structural equations of motions. The new 
locations of the interior grid points then are determined by 
solving the static equilibrium equations which result from 
a summation of forces at each grid point in the x, y and 
z coordinate directions. These static equilibrium equations 
are solved using a predictor-corrector method. The new grid 
point locations are first predicted by an extrapolation from 
the previous two time levels and then corrected using sev- 
eral Jacobi iterations of the static equilibrium equations. In 
previous Euler calculations, four Jacobi iterations were suf- 
ficient to move the mesh. 11 ' 13 For the current Euler and 
Navier-Stokes calculations, additional Jacobi iterations were 
required to move the mesh. The Euler calculations required 
up to 8 iterations, and the Navier-Stokes calculations re- 
quired up to 12 iterations. 

Because the dynamic mesh is modeled using structural 
equations, this model must represent a realistic structure. 
There is one case where this has been found to be a prob- 
lem. In a C-H-type mesh for an isolated wing, chordwise 
C-type meshes are stacked along the span to form the three- 
dimensional mesh. As these C-type meshes transition from 


the wing surface to the exterior flow field the airfoil section 
is collapsed to zero thickness to form a “tip wake” surface. 
The points around the “leading-edge” of this tip wake sur- 
face must negotiate a 360° turn. The structural equations 
modeling the mesh in this area lack stiffness in one direc- 
tion. This, in turn, causes poor convergence rates for the 
predictor-corrector procedure. It was found that specifying 
the location of the leading-edge points for the tip wake alle- 
viates this problem. When using a C-H-type mesh, another 
problem occurs for Navier-Stokes applications. In this case, 
the flow solution is very sensitive to the angle of intersection 
between the trailing edge of the airfoil and the trailing-edge- 
wake line of points. It was found that, to obtain accurate 
unsteady solutions, the slope of the trailing-edge wake must 
match the slope of the instantaneous airfoil camber line at 
the trailing edge. These slopes were matched by specifying 
the points in the wake using a quadratic function. 

AGARD Wing 445.6 


The first AGARD standard aeroelastic configuration for 
dynamic response. Wing 445. 6, 30 was tested in the 16-foot 
Transonic Dynamics Tunnel (TDT) at the NASA Langley 
Research Center. 31 The wing had a quarter-chord sweep 
angle of 45°, a panel aspect ratio of 1.65, a taper ratio of 
0.66, and a NACA 65A004 airfoil section. Several mod- 
els of the wing were tested in the TDT including full span 
and semi-span models. The model used in this study is a 
semi-span, wall-mounted model which was constructed of 
laminated mahogany. The root chord of this model was 
1.833 feet and the semi-span was 2.5 feet. In order to ob- 
tain flutter data for a wide range of Mach number and density 
conditions in the TDT, holes were drilled through the ma- 
hogany wing to reduce its stiffness. The aerodynamic shape 
of the original wing was preserved by filling these holes 
with rigid foam plastic. Flutter data for this model tested in 
air are reported in Ref. 31 over a range of Mach number 
from 0.338 to 1.141. Natural boundary-layer transition was 
allowed throughout the test. The semi-span model was at- 
tached directly to the wind tunnel wall (no splitter plate was 
used); therefore, the wing root was immersed in the wall 
boundary layer. Ref. 30 indicates that the displacement 
thickness of the wall boundary layer was 0.8 inch or less. 

As in Ref. 13, the first four natural vibration modes 
of the wing are used to model the wing structure. These 
modes, numbered 1 through 4, represent first bending, first 
torsion, second bending, and second torsion, respectively, as 
determined by a finite element analysis. 30 The modes have 
natural frequencies of 9.6 Hz for the first bending mode, 
38.17 Hz for the first torsion mode, 48.35 Hz for the second 
bending mode, and 91.54 Hz for the second torsion mode as 
determined from a ground vibration test. Because the max- 
imum amplitudes of motion were small for these tests, the 
motion due to the first four vibration modes was dominated 
by out-of-plane (or vertical) displacements. Therefore in the 
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Figure 3 Partial view of the 193 x 41 x 65 computational 
grids on the wing surface and symmetry plane. 

computational analysis, the modal deflections are defined 
only by their out-of-plane deflections. 

Results and Discussion 

Steady and unsteady Euler and Navier-Stokes results 
are computed at free-stream Mach numbers M . . = 0.96 and 
1.141. The steady-state viscous calculations are performed 
at the experimental flutter conditions which corresponded to 
a free-stream Reynolds number of 364,600 per foot for the 
Mco = 0.96 calculation and a free-stream Reynolds num- 
ber of 537,700 per foot for the M .. = 1.141 calculation. 
As mentioned in a previous section, the GAF’s at each 
free-stream Mach number also are computed using the ex- 
perimental flutter conditions. Two cases of viscous time- 
marching flutter calculations at M . = 1.141 are considered. 
In one case, the Reynolds number is specified as the exper- 
imental value for calculations at each dynamic pressure. In 
the other case, the Reynolds number is varied with the dy- 
namic pressure. 



Figure 4 Planform view of computational grids. 

The Euler and Navier-Stokes computations are per- 
formed using a 193 x 41 x 65 C-H-type grid with 193 points 
wrapped around the wing and its wake (129 points on the 
wing surface), 65 points distributed from the wing root to 
the spanwise boundary (41 points on the wing surface), and 
41 points distributed radially from the wing surface to the 
outer boundary. This mesh topology is chosen rather than 
the O-O-type topology because the wind tunnel model has a 
sheared-off tip. Both grids extend 6 root chord lengths from 
the wing to the upstream boundary, 7 root chord lengths 
from the wing to the downstream boundary, 6 root chord 
lengths from the wing to the upper and lower boundaries, 
and 1 semi-span length from the tip to the side boundary. 
A partial view of the surface mesh on the wing and sym- 
metry plane for the Euler and the Navier-Stokes grids is 
shown in Figs. 3(a) and (b), respectively. The Euler and 
Navier-Stokes grids have identical chordwise and spanwise 
distributions. Therefore, the surface mesh for both grids is 
the same. A planform view of the surface mesh for both the 
Euler grid and the Navier-Stokes grid is shown in Fig. 4. 
The distribution of points in the radial direction is not the 
same for the two grids. For the Euler grid, the grid spacing 
normal to the surface is 0.1 percent of the local chord. For 
the Navier-Stokes grid, the grid spacing normal to the sur- 
face is varied over the chord such that at least one point will 
be in the laminar sublayer of the turbulent boundary layer. 
An additional finer Navier-Stokes mesh is utilized to com- 
pute steady-state results in a grid density study. This mesh 
is a 265 x 81 x 65 C-H-type mesh identical in topology to 
the 193 x 41 x 65 mesh. The fine mesh has 177 points 
wrapped around the wing chord and 41 points distributed 
from the wing root to tip. 

For all of the calculations, the Euler and Navier-Stokes 
equations are solved using flux-vector splitting and a smooth 
flux-limiter. Convergence to steady-state is accelerated us- 
ing local time-stepping, mesh sequencing and multi-grid cy- 
cling. For time-marching calculations, the nondimensional 
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(a) Moo = 0.96. 

Figure 5 Comparison of steady-state pressure 
coefficient distributions between the 265 x 81 x 65 
C-H mesh and the 193 x 41 x 65 C-H mesh. 

global time step (based on the root chord and the free-stream 
speed of sound) is 0.03273 for both the Euler and the Navier- 
Stokes calculations. Viscous calculations are performed us- 
ing the Baldwin-Lomax turbulence model with an adiabatic 
wall temperature. Since no information is available on the 
experimental transition location, turbulence is modeled over 
the entire wing surface. This is a good approximation for 
transition locations near the leading edge of the wing. 

Steady-State Results 

Before beginning the flutter analyses of the wing, a grid 
density study was performed for the steady-state viscous 
flow fields. Results from the 265 x 81 x 65 C-H mesh 
with approximately 1.4 million grid points were used as a 
bench mark. After some investigation with different grid 
densities, the 193 x 41 x 65 C-H mesh with 0.51 million 
points was found to predict very similar steady-state surface 
pressure distributions as the finer mesh for M . = 0.96 and 
for Moo = 1-141. Figures 5(a) and (b) show the surface 
pressure coefficients at five span stations corresponding to 
the 0.0 percent, 26.0 percent, 50.4 percent, 75.5 percent. 


265x81x65 

193x41x65 




(b) Moo = 1.141. 


Figure 5 Concluded. 

and 96.0 percent of the semispan for both grids. The surface 
pressure plots indicate that the coarser grid resolves the flow 
features across the span of the wing at least as well as the 
finer grid. For this reason, the 193 x 41 x 65 grid was chosen 
for this study, and all of the subsequent results presented in 
this paper were computed using this grid to minimize the 
cost of the aeroelastic computations. 

Steady-state flowfields are used as initial conditions for 
the unsteady calculations. Comparisons of steady-state pres- 
sure coefficient contours of these initial flowfields on the up- 
per wing surface are shown in Figs. 6 and 7 to illustrate the 
basic flow characteristics of the time-marching calculations 
at Moo = 0.96 and 1.141 and a = 0°. Regions of supercrit- 
ical flow (indicated by the shaded areas) are determined in 
Figs. 6(a) and (b) by the critical pressure coefficient contour 
(C* = -0.0697). Although C* is derived from isentropic 
relations, it is used in this case to give an approximate in- 
dication of supersonic flow over the wing surface. For the 
Navier-Stokes calculation, C* is an indication of supersonic 
flow at the edge of the boundary layer based on the assump- 
tion that for an attached flow the gradient of pressure through 
the boundary layer is zero. The pressure coefficient contours, 
shown in Figs. 6(a) and (b), indicate that at = 0.96 a 
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(a) Euler. 




(b) Navier-Stokes. 


Figure 6 Comparison of steady-state pressure 
coefficient contours on the upper surface of 
Wing 445.6 at = 0.96 and a = 0°. 

large area of supercritical flow has formed on the forward 
portion of the wing. This area of supercritical flow, how- 
ever, does not terminate with a shock. Although, the Eu- 
ler calculation predicts higher levels of acceleration over the 
wing (lower minimum pressure coefficient) and a more rapid 
recompression on the inboard portion of the wing, the Eu- 
ler and Navier-Stokes calculations predict similar regions of 
supercritical flow. 

Regions of supercritical flow are determined for 
Mco = 1.141 in Figs. 7(a) and (b) by the critical pressure 


Figure 7 Comparison of steady-state pressure 
coefficient contours on the upper surface of 
Wing 445.6 at = 1.141 and a = 0°. 

coefficient contour (C* = 0.2057). At M . = 1.141, the su- 
percritical region encompasses the entire wing surface ex- 
cept for a small area around the leading edge of the wing. 
Pressure coefficient contours shown in Fig. 7(a) indicate that 
for the inviscid computation, an oblique shock has formed 
on the aft portion of the wing. On the outboard portion of 
the wing, the shock is located at approximately 75 percent 
of the local chord. Pressure coefficient contours shown in 
Fig. 7(b) indicate that for the viscous computation, a weaker 
oblique shock forms on the outboard portion of the wing at 
approximately 70 percent of the local chord. A comparison 
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Figure 8 Comparison of generalized aerodynamic forces 
for the rigid plunge and pitch of Wing 445.6 at \k. = 0.96 
and o = 0° computed using the Navier-Stokes equations. 

of Figs. 6 and 7 indicates that viscosity has a greater effect 
on the surface pressures at . = 1.141 than at XL. = 0.96. 
Significant changes in the steady-state flow conditions from 
Mco = 0.96 to 1.141 are also indicated by the formation of 
a shock at the tip of the wing. This range of free-stream 
Mach number corresponds to the range of Mach numbers 
where the experimental flutter boundary and the computed 
flutter boundary in Ref. 13 rapidly rise (See Fig. 1.) 

Pulse Transfer-Function Results 

Rigid pitching and plunging The generalized forces 
for the wing at M . = 0.96 are computed for the viscous 
case using the pulse transfer-function analysis method de- 
scribed in a preceding section. The pulse calculations are 
restarted from a steady-state flow condition at a = 0 °. 
A plunging motion and a pitching motion about the root 
quarter-chord, which are defined as modes h and 6 respec- 
tively, are analyzed. These simple modes were chosen so 


that the motion of the wing could be simulated not only by 
the deforming mesh algorithm but also by a rigid transla- 
tion and rotation of the grid. The maximum pitch amplitude 
is 1 °, the maximum plunge amplitude is 0.01 root chord 
lengths, and the frequency resolution of the pulse analysis 
is 0.1. The results of the pulse analysis for the plunging 
and pitching motions, shown in Fig. 8 , are plotted as real 
and imaginary components of the unsteady forces, . 1 h , as 
a function of the reduced frequency k based on wing root 
semichord. These generalized forces represent work divided 
by the dynamic pressure and c 2 . The generalized force Ahh 
is the lift due to plunge, A/ lH is the lift coefficient due to 
pitching, Agh is the pitching moment due to plunge, and . \ HH 
is the pitching moment due to pitch. As shown in Fig. 8 , the 
forces are independent of the way in which the mesh was 
moved which verifies the implementation and performance 
of the dynamic mesh algorithm for the Navier-Stokes com- 
putations. A similar analysis is performed in Ref. 13 to 
verify the implementation and performance of the dynamic 
mesh algorithm for Euler computations. 


Structural modes The generalized forces for the wing 
at Mjx, = 0.96 and 1.141 are computed for the inviscid and 
viscous cases using the pulse transfer-function analysis. The 
pulse calculations are restarted from a steady-state flow con- 
dition at a = 0°. The first four structural modes are ana- 
lyzed. The maximum amplitude of the generalized displace- 
ment is 0.02 for mode 1, 0.005 for mode 2, 0.02 for mode 
3, and 0.005 for mode 4. These displacements result in a 
maximum deflection of approximately 0.50 inch at the wing 
tip. The results of the pulse analyses, shown in Figs. 9 and 
10 , are plotted for the real and imaginary components of the 
unsteady forces . 1 q as functions of the reduced frequency. 
The frequency resolution of the pulse analyses is 0.04. 

A comparison of inviscid and viscous results in Fig. 
9 for XL. = 0.96 indicates slight changes in the diagonal 
terms An and An?- The remaining diagonal terms, A 33 
and A44, indicate a more significant shift in the real parts 
such that the Navier-Stokes forces are smaller in magnitude 
in comparison with the Euler forces. In general, the real 
components of the GAF’s show a greater sensitivity to 
viscosity than do the imaginary components. Figure 10 
indicates that the same is true for the GAF’s at XL. = 1.141. 
However, in contrast to the results for XL. = 0.96, Fig. 10 
indicates that the diagonal term A?? shows a more significant 
difference in the real components of the viscous and inviscid 
forces at the experimental flutter reduced frequency of 0.1. 
A decrease in the magnitude of the real component of An? 
is associated with an aft shift of the aerodynamic center 
of the wing . 16 Therefore the results shown in Fig. 10 
indicate that the effect of viscosity is to delay the aft motion 
of the aerodynamic center. This is consistent with the 
results shown in the steady-state surface pressure coefficient 
contours for XL. = 1.141. The diagonal terms An, -4 33 , 
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Figure 9 Comparison of generalized aerodynamic forces for the first 
four structural modes of Wing 445.6 at \L. = 0.96 and a = 0 °. 


and A 44 in Fig. 10 show similar effects due to viscosity as 
shown in Fig. 9 for these GAF’s. 

A comparison of results in Figs. 9 and 10 shows the 
effects of increasing the free-stream Mach number. The di- 
agonal terms An and A 33 exhibit only slight changes with 
increase in M . . However, the real component of A 22 ex- 
hibits a decrease from M . = 0.96 to 1.141 which is greater 
for the GAF’s computed from the Euler calculation. As 
noted, the decrease in the real component of A 22 indicates 
an aft shift in the aerodynamic center of the wing. The real 
component of A44 exhibits a change in sign with increasing 
\l.. . The off-diagonal term A 21 exhibits a decrease in the 
magnitude of real component with an increase in M w which 
is also an indication of an aft shift in the aerodynamic cen- 
ter. Other off-diagonal terms exhibit changes in sign with 
increase in M K . 

Aeroelastic Results 

V-g analysis Flutter characteristics are determined in 


a V-g analysis for \L. = 0.960 and 1.141 at q- = 0° using 
the harmonic loads shown in Figs. 9 and 10. Since the wing 
has a symmetric airfoil section, at zero degrees angle of 
attack, this configuration has no static aeroelastic deflections. 
Therefore, static aeroelastic characteristics did not need to 
be considered. 

The linear stability analysis was performed for three 
modal configurations: modes 1 and 2 (2 modes), modes 1 , 
2 and 3 (3 modes), and modes 1, 2, 3 and 4 (4 modes). The 
flutter characteristics determined from these V-g analyses are 
shown in Figs. 11 and 12. The computed flutter character- 
istics in terms of flutter speed index , Uf and nondimen- 

r bio a /j, 

sional flutter frequency ratio — are shown in Figs. 1 1 and 
12 along with the experimental results. An additional case 
with structural damping (</, ) was computed for the 4-mode 
analyses, but since no structural damping was measured for 
the wind-tunnel model, a estimated value of 0.03 was ap- 
plied to each mode in order to evaluate the effects of struc- 
tural damping. The results for M . = 0.960 shown in Fig. 
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Figure 10 Comparison of generalized aerodynamic forces for the first 
four structural modes of Wing 445.6 at M . = 1.141 and a = 0°. 


11 indicate that modes 3 and 4 do not significantly affect 
the flutter speed and frequency for either the Euler or the 
Navier-Stokes results. Figure 1 1 also indicates that includ- 
ing the effects of viscosity increases the flutter speed index 
by 9.6 percent for the 4-mode analysis which improves cor- 
relation with the experimental data. The addition of struc- 
tural damping stabilizes the wing which further improves 
the correlation. The flutter speed index computed from the 
Navier-Stokes GAF’s with structural damping included are 
within 5 percent of the experimental value of flutter speed 
index. The 4— mode Euler and Navier-Stokes results (with 
structural damping) are shown in Fig. 13 along with the ex- 
perimental boundary and the flutter boundary predicted in 
Ref. 13. The present 4-mode Euler flutter analysis predicts 
a lower flutter speed (approximately 7 percent) than that pre- 
dicted in Ref. 13. This is thought to be due to the increase in 
spatial resolution of the Euler grid used in the present analy- 
ses. The Reynolds numbers at flutter predicted in the V-g 
analyses with the Navier-Stokes GAF’s range from 374,000 


per foot to 382,000 per foot, a maximum of 4.5 percent dif- 
ference from the experimental free-stream Reynolds number 
(Q/Qe, P = 0.86-0.92). 

The results shown in Fig. 12 indicate that modes 3 and 
4 have a greater effect on the flutter characteristics for the 
Euler and Navier-Stokes results at \L. = 1.141. Based on 
the results of the 2-mode analysis, the significant decrease in 
the real component of . 1 2 2 noted in the Euler pulse analyses 
from M . = 0.96 to 1.141 provides the stabilizing aerody- 
namics that produced the rapid rise on the flutter boundary 
in this Mach number range. Including modes 3 and 4 in 
the analysis tends to further stabilize the wing. A compar- 
ison between the flutter characteristics predicted with the 
Euler and Navier-Stokes GAF’s in Fig. 12 indicates that in- 
cluding the effects of viscosity significantly improves the 
correlation with the experimental results. In the 4-mode 
analysis, the flutter speed index predicted with the Navier- 
Stokes GAF’s is lower than the flutter speed index predicted 
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Figure 1 1 Comparison of computed flutter characteristics 
for Wing 445.6 at XI. = 0.96 and a = 0°. 

with the Euler GAF’s by 46 percent of the experimental 
value. Figure 12 also shows that the inclusion of structural 
damping destabilizes the wing at this Mach number which 
further improves the correlation with the experiment. (The 
fact that the addition of structural damping destabilizes the 
wing is somewhat unusual). However, the flutter speed in- 
dex computed using the Navier-Stokes GAF’s with structural 
damping is 17.5 percent higher than the experimental value. 
The 4-mode Euler and Navier-Stokes results (with structural 
damping) are shown in Fig. 13 along with the experimen- 
tal boundary and the flutter boundary predicted in Ref. 13. 
The present Euler flutter analysis predicts flutter at a slightly 
higher speed (approximately 3 percent) than that predicted 
in Ref. 13. This is also thought to be due to the increase in 
spatial resolution of the Euler grid used in the present analy- 
ses. The Reynolds numbers at flutter predicted in the V-g 
analyses with the Navier-Stokes GAF’s range from 486,000 
per foot to 502,000 per foot, a maximum of 9.6 percent dif- 
ference from the experimental free-stream Reynolds number 
(Q/Qe, P = 1.28-1.46). 
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Figure 12 Comparison of computed flutter characteristics 
for Wing 445.6 at XI.. =1.141 and a = 0°. 

Time-marching analysis Time-marching calculations 
are performed at .1/ . = 1.141 using the Euler and Navier- 
Stokes equations because more significant discrepancies 
were noted between the computed and experimental results 
at this free-stream Mach number than at XI. = 0.96. In 
the time-marching analyses, the first four modes are used to 
model the structure. As mentioned previously, the viscous 
time-marching calculations are performed with and without 
variation in the free-stream Reynolds number at dynamic 
pressures ranging from Q/Q e .vp = 1-4 to 1.85 (flutter dy- 
namic pressure predicted at Q/Qe.vp = 1-61). The flutter 
results are, however, independent of this change (less than 
1 percent difference). The results of the Euler and Navier- 
Stokes time-marching aeroelastic analyses are included with 
the V-g results shown in Fig. 12. The flutter speed index 
and flutter frequency from the Euler time-marching analysis 
compare closely with the results of the Euler V-g (4-mode) 
analysis. The flutter characteristics predicted by the Navier- 
Stokes time-marching analysis are high in comparison with 
the Navier-Stokes V-g (4-mode) results by approximately 5 
percent. An additional Navier-Stokes time-marching calcu- 
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lation performed at M . = 1.141 included viscous structural 
damping added to all four modes. The magnitude of the vis- 
cous damping corresponded to the amount of gi added in the 
V-g analysis. The addition of the viscous type of structural 
damping to the time-marching analysis predicted a drop in 
the flutter speed index which is the same percentage drop as 
that predicted from the V-g analysis. 

Although the additional modeling of the viscous bound- 
ary layer did improve the correlation with the experiment for 
the cases considered in this study, several issues concerning 
the modeling of this aeroelastic problem still need to be ad- 
dressed. Since the wind tunnel test was performed under 
free-transition conditions, information concerning transition 
over the wing model is unknown, and therefore, it is difficult 
to evaluate the accuracy of the modeling of the boundary 
layer. In addition, accurate modeling of the tip geometry 
(the wind tunnel model has a cut-off tip) is difficult for 
body-fitted meshes, and the effect of this loss in geometric 
definition on the tip aerodynamics is unknown. The results 
of the V-g analysis also indicate that higher modes might be 
needed in the computational model at the higher M m . The 
effects of in-plane deflections are an additional unknown. 

Conclusions 

The flutter characteristics of an isolated 45° swept-back 
wing. Wing 445.6, were studied using an unsteady Euler 
and Navier-Stokes algorithm in order to investigate a pre- 
viously noted discrepancy between Euler flutter character- 
istics and the experimental data. The algorithm, which is 
a three-dimensional, implicit, upwind Euler/Navier-Stokes 
code (CFL3D Version 2.1), was previously modified for the 
time-marching, aeroelastic analysis of wings using the un- 
steady Euler equations. These modifications included the in- 
corporation of a deforming mesh algorithm and the addition 
of the structural equations of motion for their simultaneous 
time integration with the governing flow equations. In this 
paper, the aeroelastic method was extended and evaluated for 
applications using Navier-Stokes aerodynamics. The paper 
presented a brief description of the aeroelastic method and 
presented unsteady calculations which verified the method 
for Navier-Stokes calculations. A linear stability analysis 
and a time-marching aeroelastic analysis were used to deter- 
mine the flutter characteristics of the wing. Effects of fluid 
viscosity, structural damping and number of modes in the 
structural model were investigated. 

In the linear stability analysis, the unsteady general- 
ized aerodynamic forces of the wing were computed using 
the Euler and Navier-Stokes equations for a range of re- 
duced frequency using the pulse transfer-function method. 
The flutter characteristics of the wing were determined us- 
ing these unsteady generalized aerodynamic forces in a tra- 
ditional V-g flutter analysis at free-stream Mach numbers of 
0.96 and 1.141. The results of the V-g analysis indicated 
that the effect of viscosity on the flutter characteristics was 


Flutter 

Speed 

Index 




Mach Number 

Figure 13 Comparison of Euler and Navier-Stokes 
V-g results with Euler flutter results from Ref. 

13 and experimental data for Wing 445.6. 

to delay the rise in the flutter boundary on the supersonic 
side which significantly improved correlation with the ex- 
perimental boundary. The V-g analysis also indicated that 
the effect of structural damping was to delay the rise in the 
boundary as well. Although the flutter characteristics pre- 
dicted by the Navier-Stokes V-g analysis for the free-stream 
Mach number of 0.96 compared very well with the experi- 
mental data, the results for the free-stream Mach number of 
1.141 were still high in comparison with the experimental 
data. The effect on the flutter characteristics of including 
higher structural modes in the V-g analysis was less signif- 
icant at the free-stream Mach number of 0.96 than at the 
free-stream Mach number of 1.141. 

Time-marching aeroelastic calculations were performed 
at a free-stream Mach number of 1.141 using the Euler 
and Navier-Stokes equations to compare with the linear V- 
g flutter analysis method. The Euler results for the time- 
marching analysis were within 2 percent of the flutter speed 
predicted by the V-g analysis, and the Navier-Stokes results 
for the time-marching analysis were within 5 percent of the 
flutter speed predicted by the V-g analysis. Although the 
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linear stability analyisis flutter results did not match the time- 
marching results, the V-g analysis did, however, indicate 
the significant effect of viscosity on the supersonic flutter 
boundary for this wing. 
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